#pragma once

#include "simplify3.h"
#include <limits>

namespace outline {

template <typename E>
double simplifier<E>::simplify_min_step_len_core(int n, ///< min_step_len
	sample_vector const & src, patch_vector & dst, bool closed)
{
	assert (n >= estimator.min_step_len());

	size_t const m = src.size();
	assert (m > 1);

	here_vector ns(m);

	ns[0].d = 0;
	ns[0].p.istart = 0;

	size_t i;
	estimator.start(src[0]);
	for (i = 1; i < n && i < m; ++i) {
		ns[i].d = -1;
		ns[i].p.istart = 0;
		estimator.add_sample(src[i]);
	}
	for (; i < 2*n && i < m; ++i) {
		ns[i].d = estimator.add_sample_and_estimate(src[i], &ns[i].p);
		ns[i].p.istart = 0;
	}
	for (; i < m; ++i) {
		size_t max_j = (i < 3*n) ? i-n+1 : 2*n;
		estimator.start(src[i]);
		ns[i].d = std::numeric_limits<double>::max();
		size_t j;
		for (j = 1; j < n; ++j)
			estimator.add_sample(src[i-j]);
		for (; j < max_j; ++j) {
			assert(ns[i-j].d >= 0);
			patch p;
			p.istart = i - j;
			double d = ns[i-j].d + estimator.add_sample_and_estimate(src[i-j], &p, true);
			if (d < ns[i].d) {
				ns[i].d = d;
				ns[i].p = p;

				//test for better start
				if (closed && j == 2*n-1) {
					int s = i - j;
					int c = 1;
					for (int k = 1; k < n; ++k)
						if (ns[i-k].p.istart == s)
							++c;
						else
							break;
					if (c == n) {
						//we found better fixed start point for closed curve
						ns.clear();
						sample_vector tmp;
						shift(src, tmp, s);
						double d = simplify_min_step_len_core(n, tmp, dst, 
							false //no closed flag because second start reselection is unwanted
						);
						back_shift(dst, s, m);
						return d;
					}
				}
			}
		}
		assert(ns[i].d != std::numeric_limits<double>::max());
	}

	//back pass
	patch_vector tmp;
	for (i = m-1; i != 0; i = ns[i].p.istart) {
		tmp.push_back(ns[i].p);
	}
	dst.assign(tmp.rbegin(), tmp.rend());

	return ns[m-1].d;
}

template <typename E>
double simplifier<E>::simplify_min_step_len(int min_step_len, 
	sample_vector const & src, patch_vector & dst)
{
	if (min_step_len < estimator.min_step_len())
		min_step_len = estimator.min_step_len();

	return simplify_min_step_len_core(min_step_len, src, dst, closed(src));
}

template <typename E>
double simplifier<E>::simplify_many_min_step_len(int min_step_len,
	some_sample_vectors const & src, some_patch_vectors & dst)
{
	dst.resize(src.size());
	double s = 0;
	for (int i = 0; i < src.size(); ++i)
		s += simplify_min_step_len(min_step_len, src[i], dst[i]);
	return s;
}

//

template <typename E>
void simplifier<E>::make_first_layer(sample_vector const & src, here_vector & first)
{
	int i;
	size_t m = first.size();

	first[0].d = -1;
	estimator.start(src[0]);
	assert (estimator.min_step_len() < m);
	for (i = 1; i+1 < estimator.min_step_len(); ++i) {
		estimator.add_sample(src[i]);
		first[i].d = -1;
	}
	for (; i+1 < m; ++i) {
		first[i].d = estimator.add_sample_and_estimate(src[i], &first[i].p);
		first[i].p.istart = 0;
	}
	if (closed(src))
		first[m-1].d = estimator.add_start_and_estimate(src[m-1], &first[m-1].p);
	else
		first[m-1].d = estimator.add_sample_and_estimate(src[m-1], &first[m-1].p);
	first[m-1].p.istart = 0;
}

template <typename E>
void simplifier<E>::make_next_layer(sample_vector const & src,
	here_vector const & last, here_vector & next)
{
	int i, j;

	size_t m = last.size();
	next.resize(m);
	next[0].d = -1;
	assert (estimator.min_step_len() < m);
	for (i = 1; i+1 < estimator.min_step_len(); ++i)
		next[i].d = -1;
	for (; last[i-estimator.min_step_len()+1].d < 0 && i < m; ++i)
		next[i].d = -1;
	for (; i < m; ++i) {
		estimator.start(src[i]);
		next[i].d = std::numeric_limits<double>::max();
		for (j = i-1; i-j+1 < estimator.min_step_len(); --j)
			estimator.add_sample(src[j]);
		for (; last[j].d >= 0; --j) {
			patch p;
			p.istart = j;
			double d = last[j].d + estimator.add_sample_and_estimate(src[j], &p, true);
			if (d < next[i].d) {
				next[i].p = p;
				next[i].d = d;
			}
			else if (estimator.min_discrepancy() > next[i].d)
				break;
		}
		assert(next[i].d != std::numeric_limits<double>::max());
	}
}

///builds optimal coutour in reverse order by passing through layers
template <typename E>
void simplifier<E>::layers_back_pass(
	std::vector<here_vector> const & layers, int last_step,
	patch_vector & reverse)
{
	int i = last_step;
	int k = layers.size() - 1;
	while (i != 0) {
		reverse.push_back(layers[k][i].p);
		i = layers[k][i].p.istart;
		--k;
	}
	assert(k == -1);
}

template <typename E>
double simplifier<E>::simplify_to_K_nodes_core(int K,
	sample_vector const & src, patch_vector & dst, bool closed)
{ 
	assert (K >= 2);
	size_t const m = src.size();
	assert (m >= 2);

	int i, k;

	if (K == 2) {
		//solution is known, we have to calculate discrepancy only
		estimator.start(src[0]);
		for (i = 1; i < m-1; ++i)
			estimator.add_sample(src[i]);

		dst.resize(1);
		dst[0].istart = 0;
		return estimator.add_sample_and_estimate(src[m-1], &dst[0]);
	}

	std::vector<here_vector> layers(K-2);

	layers[0].resize(m);
	make_first_layer(src, layers[0]);

	for (k = 1; k < K-2; ++k) 
		make_next_layer(src, layers[k-1], layers[k]);

	//last step
	patch ls;
	double lsd = std::numeric_limits<double>::max();
	estimator.start(src[m-1]);
	for (i = m-2; m-i < estimator.min_step_len(); --i)
		estimator.add_sample(src[i]);
	for (; layers[K-3][i].d >= 0; --i) {
		patch p;
		p.istart = i;
		double d = layers[K-3][i].d + estimator.add_sample_and_estimate(src[i], &p, true);
		if (d < lsd) {
			ls = p;
			lsd = d;
		}
	}
	assert(lsd != std::numeric_limits<double>::max());

	if (closed && K >= 3) {
		//how many times node is leaved in contour,
		//summed for simplified contours of size 3..K
		std::vector<int> stat(m);
		for (int j = 1; j <= K-3; ++j) {
			if (j == K-3)
				i = ls.istart;
			else
				i = layers[j+1][m-1].p.istart;
			k = j;
			while (i != 0) {
				++stat[i];
				i = layers[k][i].p.istart;
				--k;
			}
			assert(k == -1);
		}
		int s = std::max_element(stat.begin(), stat.end()) - stat.begin();

		//new fixed start point for closed curve is choosen
		stat.clear();
		layers.clear();
		sample_vector tmp;
		shift(src, tmp, s);
		double d = simplify_to_K_nodes_core(K, tmp, dst, 
			false //no closed flag because second start reselection is unwanted
		);
		back_shift(dst, s, m);
		return d;
	}

	//back pass
	patch_vector tmp;
	tmp.push_back(ls);
	layers_back_pass(layers, ls.istart, tmp);
	dst.assign(tmp.rbegin(), tmp.rend());

	return lsd;
}

template <typename E>
double simplifier<E>::simplify_to_K_nodes(int K,
	sample_vector const & src, patch_vector & dst)
{
	if (K <= 1)
		K = 2;
	return simplify_to_K_nodes_core(K, src, dst, closed(src));
}

template <typename E>
double simplifier<E>::simplify_many_avg_step_len(float avg_step_len,
	some_sample_vectors const & src, some_patch_vectors & dst)
{
	dst.resize(src.size());
	double s = 0;
	for (int i = 0; i < src.size(); ++i) {
		int K = src[i].size() / avg_step_len + 0.5f;
		s += simplify_to_K_nodes(K, src[i], dst[i]);
	}
	return s;
}

//

template <typename E>
double simplifier<E>::simplify_max_discrepancy_core(double max_discrepancy,
	sample_vector const & src, patch_vector & dst, bool closed)
{
	assert(max_discrepancy > 0);
	size_t const m = src.size();
	assert (m >= 2);

	std::vector<here_vector> layers(1);
	layers.reserve(m); //no resizing is permitted

	layers[0].resize(m);
	make_first_layer(src, layers[0]);

	int k;
	for (k = 0; layers[k][m-1].d > max_discrepancy; ++k) {
		layers.push_back(here_vector());
		make_next_layer(src, layers[k], layers[k+1]);
	}

	if (closed && k >= 1) {
		//how many times node is leaved in contour,
		//summed for simplified contours of arbitrary size
		std::vector<int> stat(m);
		for (int j = 1; j <= k; ++j) {
			int i = layers[j][m-1].p.istart;
			int l = j-1;
			while (i != 0) {
				++stat[i];
				i = layers[l][i].p.istart;
				--l;
			}
			assert(l == -1);
		}
		int s = std::max_element(stat.begin(), stat.end()) - stat.begin();

		//new fixed start point for closed curve is choosen
		stat.clear();
		layers.clear();
		sample_vector tmp;
		shift(src, tmp, s);
		double d = simplify_max_discrepancy_core(max_discrepancy, tmp, dst, 
			false //no closed flag because second start reselection is unwanted
		);
		back_shift(dst, s, m);
		return d;
	}

	//back pass
	patch_vector tmp;
	layers_back_pass(layers, m-1, tmp);
	dst.assign(tmp.rbegin(), tmp.rend());

	return layers[k][m-1].d;
}

template <typename E>
double simplifier<E>::simplify_max_discrepancy(double max_discrepancy,
	sample_vector const & src, patch_vector & dst)
{
	return simplify_max_discrepancy_core(max_discrepancy, src, dst, closed(src));
}

template <typename E>
double simplifier<E>::simplify_many_avg_discrepancy(double avg_discrepancy,
	some_sample_vectors const & src, some_patch_vectors & dst)
{
	assert (avg_discrepancy >= 0);
	dst.resize(src.size());
	double s = 0;
	for (int i = 0; i < src.size(); ++i)
		s += simplify_max_discrepancy(src[i].size() * avg_discrepancy, src[i], dst[i]);
	return s;
}

template <typename E>
void simplifier<E>::back_shift(patch_vector & dst, int s, int outline_size)
{
	int pc = dst.size();
	int i;
	for (i = 0; i < pc; ++i) {
		if ((dst[i].istart += s) >= outline_size-1) {
			dst[i].istart -= outline_size-1;
			assert (dst[i].istart < outline_size);
		}
	}

	//find minimum among istart
	int j = 0;
	size_t min = dst[0].istart;
	for (i = 1; i < pc; ++i) {
		if (dst[i].istart < min) {
			j = i;
			min = dst[i].istart;
		}
	}

	patch_vector tmp(pc);
	std::copy(dst.begin() + j, dst.end(), tmp.begin());
	std::copy(dst.begin(), dst.begin() + j, tmp.begin() + pc-j);
	dst.swap(tmp);
}

} //namespace outline
